Quantitative modeling of human liver reveals dysregulation of glycosphingolipid pathways in nonalcoholic fatty liver disease

Summary Nonalcoholic fatty liver disease (NAFLD) is an increasingly prevalent disease that is associated with multiple metabolic disturbances, yet the metabolic pathways underlying its progression are poorly understood. Here, we studied metabolic pathways of the human liver across the full histological spectrum of NAFLD. We analyzed whole liver tissue transcriptomics and serum metabolomics data obtained from a large, prospectively enrolled cohort of 206 histologically characterized patients derived from the European NAFLD Registry and developed genome-scale metabolic models (GEMs) of human hepatocytes at different stages of NAFLD. We identified several metabolic signatures in the liver and blood of these patients, specifically highlighting the alteration of vitamins (A, E) and glycosphingolipids, and their link with complex glycosaminoglycans in advanced fibrosis. Furthermore, we derived GEMs and identified metabolic signatures of three common NAFLD-associated gene variants (PNPLA3, TM6SF2, and HSD17B13). The study demonstrates dysregulated liver metabolic pathways which may contribute to the progression of NAFLD.


INTRODUCTION
Nonalcoholic fatty liver disease (NAFLD) represents a histological spectrum characterized by increased lipid accumulation in the hepatocytes, that encompasses nonalcoholic fatty liver (''simple'' steatosis, NAFL) and an inflammatory form termed nonalcoholic steatohepatitis (NASH) in which progressive hepatic fibrosis occurs, ultimately leading to cirrhosis and hepatocellular carcinoma (HCC) in some patients (Anstee et al., 2019). NAFLD is associated with features of the metabolic syndrome including obesity, type 2 diabetes mellitus (T2DM), dyslipidemia, and hypertension (Anstee et al., 2013;Labenz et al., 2018). The global prevalence of NAFLD in the general population has been estimated to be 25%, while the prevalence of NASH has been estimated to range from 3% to 5% (Younossi et al., 2018). The prevalence of NAFLD has been increasing proportionately with the epidemics of obesity and is predicted to continue to rise (Estes et al., 2018).
The natural history of NAFLD is highly variable, with substantial interpatient variation in disease severity and outcome predicted by the degree of fibrosis (McPherson et al., 2015;Taylor et al., 2020). NAFLD is generally considered to be a complex disease trait driven by an obesogenic environment acting on a background of genetic susceptibility (Anstee et al., 2020). However, knowledge about the processes that contribute to the observed variation in severity remains incomplete. There is a pressing need to elucidate the pathophysiological processes that occur as NAFL progresses to NASH and fibrosis stage increases toward cirrhosis. Such knowledge will define novel biomarkers that better risk stratify patients, helping to individualize their care, and aid the development of novel therapeutic strategies and drug targets.
Progression of NAFLD is characterized by distinct metabolic changes in the liver (Masoodi et al., 2021). Previously, we showed that there is an excess accumulation of liver triacylglycerols (TGs) in NAFLD, particularly those with low carbon number and double bond content (Oresic et al., 2013), which reflects increased de novo lipogenesis in NAFLD (Kotronen et al., 2009;Westerbacka et al., 2010). iScience Article Over the past decade, genome-scale metabolic modeling (GSMM) has emerged as a powerful tool to study metabolism in human cells (Bordbar et al., 2012;Sen et al., 2017Sen et al., , 2020, including in modeling metabolism of human liver under healthy and disease conditions (Hyotylainen et al., 2016;Jerby et al., 2010;Mardinoglu et al., 2014). In a previous study, we charted metabolic activities associated with degree of steatosis in human liver from 16 NAFLD cases by integrating genome-wide transcriptomics data from human liver biopsies, and metabolic flux data measured across the human splanchnic vascular bed using GSMM (Hyotylainen et al., 2016). A similar GSMM study in 32 patients with NAFLD and 20 controls identified serine deficiency and suggested that the serum concentrations of glycosaminoglycans (GAGs) (e.g., chondroitin and heparan sulfates) associate with the staging of NAFLD . Both studies, while informative, only used transcriptomic data from relatively small patient cohorts, with crude histological characterization of disease severity, thus neither study adequately addressed grade of steatohepatitis or stage of fibrosis with sufficient granularity. In addition, several other studies in human and mice have identified molecular signatures in the liver associated with NAFLD (Arendt et al., 2015;Haas et al., 2019;Lefebvre et al., 2017;Moylan et al., 2014;Starmann et al., 2012). Moylan et al. compared the hepatic gene expression profiles in high-(severe) and low-risk (mild) NAFLD patients and idenfitied specific metabolic pathways that were differentially activated in these groups (Moylan et al., 2014). In a cross-sectional study, Arendt et al. compared the hepatic gene expression in subjects with either healthy liver, NAFL or NASH and showed marked alteration in the gene expression profiles associated with the polyunsaturated fatty acid metabolism (Arendt et al., 2015). Integrative analyisis of human transcriptomic datasets and the mice models have identified functional role of dermatopontin in collagen deposition and hepatic fibrosis (Lefebvre et al., 2017). Futhermore, transcriptomic analysis and immune profiling of NASH patients also relvealed that several genes associated with the inflammatory responses were altered in progression to NASH (Haas et al., 2019). However, these studies have largely employed microarray-based techniques and hence lacked the comprehensive approach provided by the genome-wide RNA sequencing (RNA-seq) .
Herein, we examined the changes in hepatic metabolic processes and pathways that occur across the full histological spectrum of NAFLD, exploring the evolving changes during the progression of NAFL to NASH with increasing fibrosis stages including cirrhosis (F0 up to F4). Whole liver tissue RNA-Seq transcriptomic data from a cohort of histologically characterized patients derived from the European NAFLD Registry (n = 206) Hardy et al., 2020) were used to develop personalized and NAFLD stagespecific genome-scale metabolic models (GEMs) of human hepatocytes (Tables 1, S1, and S2; Figure 1). The integrative approach employed in this study defined the changes of liver metabolism in NAFL and with progressive NASH-associated fibrosis. Furthermore, GSMM predicted the metabolic differences among carriers of widely validated genetic variants associated with NAFLD/NASH disease severity in three genes (PNPLA3, TM6SF2, and HSD17B13) (Abul-Husn et al., 2018;Anstee et al., 2020;Liu et al., 2014;Romeo et al., 2008).

Development of genome-scale metabolic models of human hepatocytes in the patients with NAFLD
A GEM of human hepatocyte (iHepatocytes2322) was first developed by Mardinoglu et al., by combining the iHepatocyte1154 GEM (Agren et al., 2012) with previously published liver models. iHepatocytes2322 comprises 2,322 genes, 5,686 metabolites, and 7,930 reactions found in the human liver . Here, by using NAFLD stage-specific transcriptomics data generated from the liver biopsies of 206 patients with NAFLD (Tables 1 and S1) , we contextualized the iHepatocytes2322 model and developed personalized GEMs of human hepatocytes across the full iScience Article spectrum of NAFLD severity (NAFL through NASH F0-F4) ( Figure 1). A step-wise strategy that combines iMAT  and E-flux (Colijn et al., 2009) algorithm was developed for the contextualization of GEMs, i.e., selection of stage-specific active reactions and their associated genes and metabolites from the iHepatocytes2322. The number of genes, metabolites, and reactions included in these models is given in ( Figure S1). All the models were tested to carry out 256 metabolic tasks as given in the study by . The occurrence scores of the metabolic tasks carried out by these models are shown in ( Figure S2).
Stage-specific alterations of the metabolic pathways at various stages of NAFLD Partial least squares-discriminant analysis (PLS-DA) (Le Cao et al., 2011) of the metabolic fluxes predicted by the personalized hepatocyte GEMs showed that patients with NASH (F3, AUC = 0.70) and cirrhosis  2C). Here, the reaction flux states were estimated by a non-uniform random sampling method that finds solutions among the feasible flux distributions of the metabolic network (Bordel et al., 2010;Osterlund et al., 2013). The PLS-DA model has identified several metabolic subsystems/pathways such as fatty acid, glycerophspholipid, cholesterol, bile acid, tryptophan, histidine, glutathione, vitamin (B,D), nucleotide metabolism, and ROS detoxification which were altered (Variable Importance in Projection (VIP) scores (Farré s et al., 2015) >1) among these groups ( Figure S3). Pairwise comparisons between different NAFLD groups suggested that most of these subsystems were upregulated in the patients with ''clinically significant'' non-cirrhotic or ''advanced'' fibrosis as compared to ''mild'' or ''minimal'' diseases ( Figure 2B). Intriguingly, at ''advanced'' (vs. ''mild'') fibrosis, fluxes across the sphingolipid metabolism were significantly (two sample t-test, p values < 0.05, adjusted for False Discovery Rate (FDR)) upregulated, while glycosphingolipid (GSL) biosynthesis-globo series was downregulated ( Figure 2B).
Reporter metabolite analysis of human hepatocytes unveiled the intrinsic regulation of GSLs, GAGs, and cholesterols in advance fibrosis Reporter metabolite (RM) analysis is a metabolite-centric differential approach that aids in identifying metabolites in a network around which significant transcriptional changes occur ( iScience Article Patil and Nielsen, 2005). RMs can predict hot spots in a metabolic network that are significantly altered between two different conditions, in this case, different stages of NAFLD.
Taken together, in ''advanced'' fibrosis, the hepatocytes tend to decreases the synthesis of GSLs, particularly GlcCers, LacCers, and GalCers from Cers ( Figures 4A-4C). Excessive accumulation of Cers in the hepatocyte may be required to exert cellular signaling, modulate inflammatory response, or generate apoptotic stress signals (Apostolidis et al., 2016;Zhang et al., 2019). To validate the RMs and pathway(s) predicted by GSMM, we measured the serum levels of Cers, GSLs, and retinoids in 360 individuals participated in the EPoS cohort ( Figure 5). 41 individuals overlapped with the present study ( Figure S10).
Overall, the directional changes (up-or downregulated) of the predicted RMs of GSLs in the liver were comparable with the differences (increase or decrease) in the serum GSL levels of the patients at different stages of NASHassociated fibrosis. Intriguingly, there was an abrupt decrease in the levels of Cers and HexCers at NASH F4 (cirrhosis) stage, which may be attributed to the alterations in the SL pathway(s) in the liver (Figures 5, S10, and S11).
We developed personalized GEMs for hepatocytes in the individuals carrying one of these three common gene variants, PNPLA3 (GC, GG) (n = 69), TM6SF2 (CT, TT) (n = 13), and HSD17B13 (-T, TT) (n = 21) compared with wild type (WTs) (n = 36), i.e., individuals homozygous for all the three gene variants: PNPLA3 (CC), TM6SF2 (CC), and HSD17B13 (-)( Table 1 and Figure S14), using genome-wide transcriptomics data and the iHepatocytes2322 model. Selection of individuals exclusively carrying PNPLA3, TM6SF2, or HSD17B13 gene variants from the study cohort  is given in (Table S3). Assignment of gene variants to the disease groups is given in (Table S4). The occurrence scores of the metabolic tasks carried out by the gene variants and WT models are shown in ( Figure S15).
Pairwise PLS-DA analysis of the fluxes across the metabolic subsystems predicted by the hepatocyte-GEMs of PNPLA3, TM6SF2, or HSD17B13 variant carriers vs. WT showed remarkable separation between these groups, with AUCs = [0.86, 0.85, and 0.84], respectively ( Figures 6A-6C). PLS-DA identified five metabolic subsystems, fatty acid biosynthesis and oxidation, oxidative phosphorylation, terpinoid biosynthesis, and folate metabolism, which were commonly altered (VIP >0.1) ( Figures 6D and S16) between the gene variants and the WT. In addition, a pairwise comparison between the PNPLA3 variants and WT groups suggests that GSL metabolism was downregulated (two-sample t-test, p < 0.05 adjusted for FDR) while cysteine and methionine, and vitamin B 6 metabolism was upregulated (p < 0.05 adjusted for FDR). PNPLA3 and HSD17B13 variants showed a decrease (p < 0.05 adjusted for FDR) in the b-oxidation of fatty acids as compared to the WT, while TM6SF2 variant showed an elevation in the cholesterol biosynthesis ( Figure 6E). However, GSMM could not predict changes in the retinol (Vitamin A) metabolism in the individuals carrying PNPLA3 gene variants, as previously reported in the study by (Kovarova et al., 2015;Mondul et al., 2015).
In order to validate the predictions from GSMM, we identified the serum metabolic profiles of the individuals exclusively carrying PNPLA3 (n = 14), TM6SF2 (n = 3), or HSD17B13 (n = 2) gene variants, and WT groups (n = 11) from a subgroup of 41 individuals/samples on which RNA-Seq and metabolomic measurement was performed. A decrease (p < 0.05 adjusted for FDR) in the levels of Cers, GSLs, and retinyl palmitate in the individuals exclusively carrying PNPLA3 (vs. WT group) was observed ( Figures 6F, 6G, and S17).

DISCUSSION
By performing GSMM of human hepatocytes in NASH-associated fibrosis, we were able to identify several metabolic signatures and pathways markedly regulated at different stages of fibrosis. The previous GSMM study of human NAFLD  also identified altered metabolic pathways in NAFLD, however, that study was confined to a relatively small patients group, not representing the full spectrum of NAFLD. Furthermore, the study was compounded by lack of lipidomic/metabolomic measurement in the patient samples.
Here, we used genome-wide transcriptomics (RNA-seq)  data that provided a dynamic range of genes, enzymes/reactions, metabolites and their interactions, to model stage-specific GEMs for human hepatocytes in NAFLD (Figure 1). In addition, the study cohort allowed us to stratify and model the metabolic differences in three major genetic variants (PNPLA3, TM6SF2, and HSD17B13) associated with risk and severity of NAFLD.
Our results implicate changes in the levels of vitamins (A, E) in the NASH-associated fibrosis progression. In the hepatocytes, a significant increase in the RMs of retinoic acid derivatives was observed in NASH F2, F3 (vs. NAFL) and in ''clinically significant'' non-cirrhotic fibrosis (vs. ''minimal disease''). A small increase in the serum levels of vitamin A (retinol) and retinyl palmitate was observed at NASH F0-1 stage, while retinyl palmitate was markedly decreased at NASH F4 stage (Figures S12 and S13 iScience Article deficiency (Pettinelli et al., 2018;Saeed et al., 2017). Vitamin A is a regulator of glucose and lipid metabolism in the liver and adipose tissue and may attenuate in the development of NAFLD (Pettinelli et al., 2018;Saeed et al., 2017). The PNPLA3 gene product is known to have retinyl ester hydrolase activity, and PNPLA3-I148M is associated with low serum retinol level with enhanced retinyl esters in the liver of patients with NAFLD (Kovarova et al., 2015;Mondul et al., 2015). However, GSMM has not identified any differences in vitamin A metabolism related to PNPLA3 variant carriage. Intriguingly, vitamin E derivatives were subsequently decreased across the stages of fibrosis. High dose of vitamin E supplementation has been shown to improve histological steatohepatitis over placebo in the PIVENS randomized controlled trial of pioglitazone or vitamin E in patients with NASH (Sanyal et al., 2010).
We observed a dynamic regulation of complex GSLs in the advanced fibrosis. At that stage, the majority of cerebrosides (HexCers: glucosylceramides (GlcCers) and lactosylceramides (LacCers)) and globosides were decreased in the patient groups. Intriguingly, the serum concentrations of GSLs (HexCers, GlcCers) were markedly decreased in the patients with advanced fibrosis or cirrhosis (NASH F4). Cers are the key intermediates of sphingolipid metabolism that promote cellular proliferation, differentiation, iScience Article and cell death (Gault et al., 2010;Turpin-Nolan and Bruning, 2020). Cers interact with several pathways involved in insulin resistance, oxidative stress, inflammation, and apoptosis that are linked to NAFLD (Gault et al., 2010;Pagadala et al., 2012). Understanding the role of Cers in the staging of NASH-associated fibrosis is of great diagnostic interest (Pagadala et al., 2012). In normal physiological conditions, Cers converts to GSLs, which prevents its excessive accumulation in the cells (Ishay et al., 2020). Our data suggest that in patients with advanced fibrosis, the serum concentrations of Cers increase up to NASH F3, and decrease at NASH F4 stages.
Differential flux analysis showed an increase of Cer production via de novo pathway, while production of GSLs (GlcCer, GalCer, and LacCers) decreased in patients with advanced fibrosis. Thus, the conversion of Cers to GSLs (via glucosylceramide synthase, GCS) might have been compromised in ''advanced'' fibrosis. Intriguingly, we found that the GSLs (LacCers and globosides) were associated with GAGs (heparin, keratin sulfate), which was also altered in the progression of NASH-associated fibrosis. It is well demonstrated that inhibiting GSL synthesis in obese mice has improved glucose homeostasis and markedly reduced the development of NAFL (Zhao et al., 2009).
Taken together, we identified several metabolic signatures and pathways in progressive stages of NAFLD, which enabled us to examine how different metabolites and their intermediates are regulated across various stages of NASH-associated fibrosis. While some of our key results corroborate with the earlier findings, others allowed us to further elucidate the dys(regulation) of metabolic pathways in NAFLD. Moreover, GSMM has demonstrated several stage-specific metabolic changes, which might help discover biomarkers, identify drug targets, and ultimately lead to effective therapeutic strategies for NASH. To this end, our data indicate the significance of GSL pathways and targeting these pathways might ameliorate the liver pathology associated with NAFL and NASH-associated fibrosis. However, some of these findings remain to be validated by in vivo and/or in vitro experiments.

Limitations of the study
We acknowledge some limitations of our study, such as the relatively small sample size, lack of metabolomics/lipidomics data from the liver biopsy samples, and small overlap of patient-matched samples from transcriptomics and metabolomics studies. Nevertheless, we report sphingolipid changes along the stepwise progression of NASH. It is, however, clear that these findings remain to be investigated in larger studies.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Matej Ore si c (matej.oresic@oru.se).

Materials availability
This study did not generate new unique reagents.

Data and code availability
Scripts for GSMM, contextualization and data analysis can be downloaded from (https://github.com/ parthoBTK/Personalized_Liver_Models.git). The personalized GEMs (.mat files) of human-hepatocytes contextualized for different stages of NAFLD, and three major gene variants are available upon request. In addition, the scripts can be downloaded from Zenodo.org: https://doi.org/10.5281/zenodo.6948647.
The RNA-Seq data are available in the NCBI GEO: GSE135251.
The lipidomic/metabolomic datasets generated in this study were submitted to the Metabolomics Workbench https://doi.org/10.21228/M85976; assigned a project ID project ID PR001095 and study ID ST001964

EXPERIMENTAL MODEL AND SUBJECT DETAILS
A full description of the transcriptomic analysis of the whole liver RNA-Seq data used in this study has previously been reported . All of the patients are from European descent. Patient samples were derived from the European NAFLD Registry (NCT04442334) (Hardy et al., 2020) and comprised snap-frozen biopsy samples and associated clinical data from 206 patients diagnosed with histologically characterized NAFLD in France, Germany, Italy, or the United Kingdom. All biopsies samples obtained in this study were centrally scored by two expert liver pathologists according to the semiquantitative NASH-Clinical Research Network 'NAFLD Activity Score' (NAS) (Kleiner et al., 2005). Fibrosis was staged from F0 to F4 (cirrhosis) . Patients with alternate diagnoses and etiologies were excluded, including excessive alcohol intake (30 g per day for males and 20 g per day for females), viral hepatitis, autoimmune liver diseases, and steatogenic medication use. As previously described (Hardy et al., 2020), collection and use of samples and clinical data for this study were approved by the relevant local and/or national Ethical Review Committee covering each participating center, with all patients providing informed consent for participation. All participant recruitment and informed consent processes at recruitment centers were conducted in compliance with nationally accepted practice in the respective territory and in accordance with the World Medical Association Declaration of Helsinki 2018.

METHOD DETAILS Transcriptomics
Transcriptomic (RNA-seq) data associated with this study was obtained from ; also available in the NCBI GEO repository: GSE135251. Differentially expressed genes between a paired conditions were either obtained , or estimated by a method stated in . All genes that were differentially expressed between the case vs. control at (p-values adjusted for FDR < 0.05) were included in the RM analysis. Moreover, the transcriptomic datasets were corrected for batch effect, gender and other confounding factors as stated in . iScience 25, 104949, September 16, 2022 iScience Article Genotyping Genotypes were obtained from the GWAS data described in detail in (Anstee et al., 2020). Some data was also obtained from RNA-seq as reported elsewhere .
The samples were analyzed by ultra-high-performance liquid chromatography quadrupole time-of-flight mass spectrometry (UHPLC-QTOFMS). Briefly, the UHPLC system used in this work was a 1290 Infinity II system from Agilent Technologies (Santa Clara, CA, USA). The system was equipped with a multi sampler (maintained at 10 C), a quaternary solvent manager and a column thermostat (maintained at 50 C). Injection volume was 1 mL and the separations were performed on an ACQUITY UPLCâ BEH C18 column (2.1 mm 3 100 mm, particle size 1.7 mm) by Waters (Milford, MA, USA). The mass spectrometer coupled to the UHPLC was a 6545 QTOF from Agilent Technologies interfaced with a dual jet stream electrospray (Ddual ESI) ion source. All analyses were performed in positive ion mode and MassHunter B.06.01 (Agilent Technologies) was used for all data acquisition. Quality control was performed throughout the dataset by including blanks, pure standard samples, extracted standard samples and control serum samples, including in-house serum and a pooled QC with an aliquot of each sample was collected and pooled and used as quality control sample.
Relative standard deviations (% RSDs) for identified lipids in the control serum samples (n = 13) and in the pooled serum samples (n = 54) were on average 22.4% and 17.5%, respectively.
Mass spectrometry data processing was performed using the open source software package MZmine 2.18 (Pluskal et al., 2010). The following steps were applied in this processing: (i) Crop filtering with a m/z range of 350-1200 m/z and an RT range of 2.0 to 12 minutes, (ii) Mass detection with a noise level of 750, (iii) Chromatogram builder with a minimum time span of 0.08 min, minimum height of 1000 and a m/z tolerance of 0.006 m/z or 10.0 ppm, (iv) Chromatogram deconvolution using the local minimum search algorithm with a 70% chromatographic threshold, 0.05 min minimum RT range, 5% minimum relative height, 1200 minimum absolute height, a minimum ration of peak top/edge of 1.2 and a peak duration range of 0.08-5.0, (v), Isotopic peak grouper with a m/z tolerance of 5.0 ppm, RT tolerance of 0.05 min, maximum charge of 2 and Stage-specific functional GEMs of human hepatocytes were developed by step-wise combining iMAT  and E-Flux (Colijn et al., 2009) algorithms, applied to iHepatocytes2322 as a template model (GEM). iMAT approach finds the optimal trade-off between inclusion and exclusion of high and low-expression reactions. It does not explicitly depends on an metabolic objective function Zur et al., 2010). iMAT requires three sets of input reactions such as high, low and moderately expressed. The reaction weights were determined by integrating the NAFLD stage-specific gene expression data to iHepatocytes2322 using gene-protein-reaction associations (GPR) rules. To categorize the model reactions into high, low and moderately expressed, we estimated the (mean G sd) of the log-normal distribution of ll OPEN ACCESS